Supplementary 3: Modelling TOC concentration and forecast in coastal drainage basins
Camille M. Crapart
24 06 2022
library(dplyr)
library(tibble)
library(DBI)
library(sf)
library(gimms)
library(ncdf4)
library(rgdal)
library(stringr)
library(ggplot2)
library(RColorBrewer)
library(colorspace)
library(cowplot)
library(png)
library(spdep)
library(spatialreg)
library(kableExtra)
library(reshape2)
library(grid)
library(gridExtra)
library(ggpubr)
source("f_plotspatial.R")
memory.limit(size = 160000)
## [1] 160000
library(captioner)
fig_nums <- captioner()
table_nums <- captioner(prefix = "Table")1 Coastal drainage basins data
1.1 Connections to database
The coastal drainage basins polygons were computed with a seamless digital elevation model for Fennoscandia (“Environmental data NINAnor/NOFA Wiki” n.d.). The model resulted in more than 200 000 coastal drainage basins, most of them being extremely small polygons covering the coastal regions. In order to simplify the model and reduce the amount of data, we keep only the 0,001 % of the coastal drainage basins that are the biggest. They account to more than 93% of the total surface covered by the initial amount of polygons.
con <- dbConnect(RPostgreSQL::PostgreSQL(),user = "camille.crapart", password = "camille",host = "vm-srv-wallace.vm.ntnu.no", dbname = "nofa")
query <- paste("SELECT gid,geom FROM ","\"Hydrography\"",".waterregions_dem_10m_nosefi",sep = "")
query.2 <- paste("SELECT * FROM ","\"Hydrography\"",".Fenoscandia_LakeRegister",sep = "")
water.regions <- st_read(con,query = query)
lake.register <- st_read(con,query = query.2)
saveRDS(water.regions,"water.regions.Rdata")water.regions <- readRDS("water.regions.Rdata")
ggplot()+geom_sf(data=water.regions,aes(fill = gid))water.regions <- readRDS("water.regions.Rdata")
water.regions$area <- st_area(water.regions) %>% as.numeric()
area1 <- sum(water.regions$area)
wr.quantile <- quantile(water.regions$area, probs = c(0.99,0.995,0.999))
wr.sf <- water.regions %>% filter(area > wr.quantile[3])
area2 <- sum(wr.sf$area)
ratio <- area2/area1*100
map.wr <- ggplot(wr.sf) + geom_sf(aes(fill=gid,col=gid))
ggsave(plot = map.wr, filename = "WR/wr.sf.png",dpi = "retina", width = 10, height = 15, units = "cm")
saveRDS(wr.sf,"WR/wr.sf.rds")
st_write(wr.sf, "wr.sf.shp")1.2 NDVI
NDVI values are extracted from the GIMMS NDVI3g dataset (The National Center for Atmospheric Research 2018), stored on http://poles.tpdc.ac.cn/en/data/9775f2b4-7370-4e5e-a537-3482c9a83d88/, in the “ecocast” dataset and accessed via the “gimms” package (Detsch 2021) on Rstudio.
Data is taken bi-monthly (Tucker et al. 2005). Raster layers of all slices are downloaded using the gimmsdownload function, then monthly composite are calculated using the monthlyComposites function. The maximum NDVI value is taken for each pixel. Afterwards the NDVI values for each catchment polygon are extracted using the extract function from the raster package (Hijmans et al. 2018).
The values stored in the ecocast dataset are composed of the NDVI value and a flag value indicating the goodness of the data (Detsch 2021). All the NDVI values extracted, corresponding to the values for the studied catchments, had a flag of 1, indicating a good value. The NDVI value was retrieved from the stored value using the formula \(NDVI = floor(ndvi3g/10)/1000\) (Detsch 2021), and then the mean of the 3 summer values (June, July and August) was computed for each catchment.
ndvi.max <- readRDS("ndvi.max.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
summer.mean <- raster::stack(ndvi.max[[6]],ndvi.max[[7]],ndvi.max[[8]]) %>% mean()
summer.ndvi <- reclassify(summer.mean,cbind(-Inf, 0, NA), right=FALSE)
summer.scandinavia <- raster::crop(summer.ndvi,c(0,35,55,73))
writeRaster(summer.mean,"ndvi1994.tif",overwrite = T)
ndvi.wr <- raster::extract(summer.ndvi,wr.sf,fun = mean, df = T, sp = T, na.rm = T)
ndvi.wr$area <- NA
ndvi.wr$ndvi <- (floor(ndvi.wr$layer)/10)/1000
saveRDS(ndvi.wr, "WR/ndvi.wr.rds")
knitr::include_graphics("WR/wr.map.NDVI.png")
Figure 1: Summer NDVI in 1994 in coastal drainage basins
1.3 Runoff
The runoff data was downloaded from the CORDEX database: https://portal.enes.org/data/enes-model-data/cordex. The data is available at https://esg-dn1.nsc.liu.se/projects/esgf-liu/.
The names of the variables in the CORDEX project follow the CF Metadata convention: http://cfconventions.org/
THe mean of the 30 years between 1970 and 2000 is used (p11 on EURO-CORDEX guidelines https://www.euro-cordex.net/imperia/md/content/csc/cordex/guidance_for_euro-cordex_climate_projections_data_use__2021-02_1_.pdf). This time period matches the temperature and precipitation data provided by WorldClim.
runoff.raster <- "CORDEX/mrros_Lmon_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_185001-201412.nc"
runoff.nc <- nc_open(runoff.raster)
runoff.stack <- raster::stack(runoff.raster, bands = c(1441:1812))
runoff.mean <- runoff.stack %>% mean()
wr.sf <- readRDS("WR/wr.sf.rds")
runoff.wr95 <- raster::extract(runoff.mean, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(runoff.wr95) <- c("gid","area","runoff")
saveRDS(runoff.wr95,"WR/runoff.wr95.rds")knitr::include_graphics("WR/wr.map.Runoff.png")
Figure 2: Runoff in 1995 in coastal drainage basins (mm/y)
1.4 Land Cover
The Land Cover data was downloaded from https://land.copernicus.eu/pan-european/corine-land-cover/lcc-2000-2006. The previous version (1995-2000) excludes Norway so the 2000 version was preferred.
The land cover data was downloaded as a raster (.tiff file), which included the legend recording the code for each land use. The raster and legend were assembled in QGIS before being imported as a shapefile in R. The codes (corresponding to the line in the extracted dataframe) of the categories of interest were:
Arable land:
- 12: non-irrigated arable land
- 16: fruit trees and berries plantations
- 18: pastures
- 19: annual crops associated with permanent crops
- 20: Complex cultivation patterns
Bogs:
- 36: peat bogs
- 35: inland marshes, excluded because all zeros in the studied area.
Forest:
- 23: Broad-leaved forest
- 24: Coniferous forest
- 25: Mixed forest
Bare:
- 31: Bare rocks
- 32: Sparsely vegetated areas
- 33: Burnt areas
wr.sf <- readRDS("WR/wr.sf.rds")
corine <- raster::raster("Corine_Land_Cover/U2006_CLC2000_V2020_20u1.tif")
wr.corine <- st_transform(wr.sf, crs(corine))
extract_corine <- function(x){
clc <- extract(corine,x)
saveRDS(clc,file = paste("WR/clc/clc",x$gid,"rds",sep="."))
}
wr.list <- split(wr.corine,seq(nrow(wr.corine)))
lapply(wr.list,extract_corine)
#clc.test <- readRDS("WR/clc/clc.16711649.rds")
clc.files <- list.files(path = "WR/clc/", pattern = "clc.*.rds", full.names = T)
list.clc <- sapply(clc.files,readRDS)
list.clc.wr <- sapply(list.clc, as.integer)
names.list <- str_extract(clc.files,"\\d+")
names(list.clc.wr) <- names.list
saveRDS(list.clc.wr,"WR/clc/list.clc.wr.rds")list.clc.wr <- readRDS("WR/clc/list.clc.wr.rds")
clc.tab.wr <- sapply(list.clc.wr, function(x) tabulate(x,45)) %>% t()
clc.tab.area.wr <- clc.tab.wr*prod(res(corine))
catchment.area <- rowSums(clc.tab.area.wr)
clc.tab.prop.wr <- sweep(clc.tab.area.wr,1,catchment.area, FUN = "/")
saveRDS(clc.tab.prop.wr,"WR/clc/clc.tab.prop.wr.rds")clc.tab.prop.wr <- readRDS("WR/clc/clc.tab.prop.wr.rds") %>% as.data.frame()
bogs <- clc.tab.prop.wr[,36] %>% as.data.frame() %>% setNames("bogs") %>% tibble::rownames_to_column(var = "gid")
saveRDS(bogs,"WR/clc/bogs.wr.rds")
arable <- rowSums(clc.tab.prop.wr[,c(12,16,18,19,20)]) %>% as.data.frame() %>% setNames("arable") %>% tibble::rownames_to_column(var = "gid")
saveRDS(arable,"WR/clc/arable.wr.rds")
forest <- rowSums(clc.tab.prop.wr[,c(23,24,25)]) %>% as.data.frame() %>% setNames("forest") %>% tibble::rownames_to_column(var = "gid")
saveRDS(forest,"WR/clc/forest.wr.rds")
bare <- rowSums(clc.tab.prop.wr[,c(31,32,33)]) %>% as.data.frame() %>% setNames("bare") %>% tibble::rownames_to_column(var = "gid")
saveRDS(bare,"WR/clc/bare.wr.rds")list.grob <- c("WR/wr.map.Bog.png", "WR/wr.map.Arable.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow = 1, ncol = 2, labels = c("a)","b)"), label_size = 10)
Figure 3: Proportion of a) bogs and b) arable land in 2000 in coastal drainage basins
1.5 Temperature and precipitation
Temperature and precipitation are downloaded from the WorldClim database: https://worldclim.org/data/index.html
Temperature in historical Worldclim is in °C * 10 (to reduce file size). Precipitation are average monthly precipitation in mm. Here we use the mean monthly precipitation and temperature for the period 1970-2000.
wr.sf <- readRDS("WR/wr.sf.rds")
bio.fennoscandia <- raster::getData(name = "worldclim", var = "bio", res = 2.5)
temp.wr <- raster::extract(bio.fennoscandia$bio1, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(temp.wr) <- c("gid","area","temp")
temp.wr$temp <- temp.wr$temp / 10
saveRDS(temp.wr,"WR/temp.wr.rds")
#temp.wr <- readRDS("WR/temp.wr.rds")
prec.wr <- raster::extract(bio.fennoscandia$bio12, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(prec.wr) <- c("gid","area","prec")
saveRDS(prec.wr,"WR/prec.wr.rds")
#prec.wr <- readRDS("WR/prec.wr.rds")
tp <- cbind(wr.sf[,1],temp.wr[,3]) %>% cbind(prec.wr[,3])
saveRDS(tp,"WR/tp.rds")list.grob <- c("WR/wr.map.Temp.png", "WR/wr.map.Precip.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow= 1, ncol = 2, labels = c("a)", "b)"), label_size = 10)
Figure 4: Map of a) average temperature (°C) and b) average precipitation (mm) in 1995 in coastal drainage basins
1.6 Nitrogen deposition (TNdep)
The nitrogen deposition data was extracted from the EMEP database (https://emep.int/mscw/mscw_moddata.html#Comp). The data from 2000 was used as it is the first model available.
N-deposition (NOx and NH3) is the sum of:
- Dry deposition of oxidized nitrogen per m2 grid DDEP_OXN_m2Grid, in mg/m2
- Wet deposition of oxidized nitrogen WDEP_OXN, in mg/m2
- Dry deposition of oxidized nitrogen per m2 grid DDEP_RDN_m2Grid, in mg/m2
- Wet deposition of reduced nitrogen WDEP_RDN, in mg/m2
library(ncdf4)
EMEP_file <- "EMEP/EMEP01_rv4.42_year.2000met_2000emis_rep2021.nc"
wr.sf <- readRDS("WR/wr.sf.rds")
woxn <- raster(EMEP_file, varname = "WDEP_OXN")
doxn <- raster(EMEP_file, varname = "DDEP_OXN_m2Grid")
wrdn <- raster(EMEP_file, varname = "WDEP_RDN")
drdn <- raster(EMEP_file, varname = "DDEP_RDN_m2Grid")
woxn_wr <- extract(woxn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(woxn_wr) <- c("gid","area","woxn")
saveRDS(woxn_wr,"woxn_wr.rds")
doxn_wr <- extract(doxn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(doxn_wr) <- c("gid","area","doxn")
saveRDS(doxn_wr,"doxn_wr.rds")
wrdn_wr <- extract(wrdn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(wrdn_wr) <- c("gid","area","wrdn")
saveRDS(wrdn_wr,"wrdn_wr.rds")
drdn_wr <- extract(drdn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(drdn_wr) <- c("gid","area","drdn")
saveRDS(drdn_wr,"drdn_wr.rds")
ndep.wr <- merge(woxn_wr,doxn_wr, by = "gid") %>% merge(wrdn_wr, by = "gid") %>% merge(drdn_wr, by = "gid")
ndep.wr$tndep <- rowSums(ndep.wr@data[,c(3:6)])
saveRDS(ndep.wr,"WR/ndep.wr.rds")knitr::include_graphics("WR/wr.map.TNdep.png")
Figure 5: Total nitrogen deposition in 1995 in coastal drainage basins (mg/m2)
1.7 Elevation
The altitude data was accessed through the “getData” function, and comes from the NASA digital elevation model (SRTM 90 m). The slope was then computed using the “terrain” function from the raster package.
wr.sf <- readRDS("WR/wr.sf.rds")
alt.nor <- raster::getData("alt",country = "NOR", mask = T)
alt.swe <- raster::getData("alt",country = "SWE", mask = T)
alt.fin <- raster::getData("alt",country = "FIN", mask = T)
alt.list <- list(alt.nor,alt.swe,alt.fin)
alt <- do.call(raster::merge,alt.list)
alt.wr <- raster::extract(alt, dplyr::select(wr.sf, c("gid","geom")), sp = T, fun = mean, df = T, na.rm = T)
names(alt.wr) <- c("gid","alt")
saveRDS(alt.wr,"WR/alt.wr.rds")
slope.nor <- terrain(alt.nor, unit = "degrees")
slope.swe <- terrain(alt.swe, unit = "degrees")
slope.fin <- terrain(alt.fin, unit = "degrees")
slope.list <- list(slope.nor,slope.swe,slope.fin)
slope <- do.call(raster::merge,slope.list)
slope.nona <- reclassify(slope,cbind(NA,100), right = NA)
slope.wr <- raster::extract(slope,wr.sf, sp = T, df = T, fun = mean, na.rm = T)
names(slope.wr) <- c("gid","area","slope")
saveRDS(slope.wr, "WR/slope.wr.rds")1.8 Gathers dataset of coastal drainage basins
The data was gathered and merged into a single dataframe, based on the catchment identifyer (“gid”).
ndvi.wr <- readRDS("WR/ndvi.wr.rds")
ndvi.wr$area <- NULL
runoff.wr95 <- readRDS("WR/runoff.wr95.rds")
runoff.wr95$area <- NULL
bogs.wr <- readRDS("WR/clc/bogs.wr.rds")
arable.wr <- readRDS("WR/clc/arable.wr.rds")
#clc.wr.sum <- readRDS("WR/clc/clc.wr.sum.Rdata")
alt.wr <- readRDS("WR/alt.wr.rds")
slope.wr <- readRDS("WR/slope.wr.rds")
slope.wr$area <- NULL
tp <- readRDS("WR/tp.rds")
ndep.wr <- readRDS("WR/ndep.wr.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
all.wr.sf <- merge(wr.sf,runoff.wr95,by = "gid") %>% merge(ndvi.wr,by="gid") %>% merge(bogs.wr,by = "gid") %>% merge(arable.wr, by = "gid") %>% merge(st_drop_geometry(tp),by = "gid") %>% merge(alt.wr, by = "gid") %>% merge(slope.wr, by = "gid") %>% merge(ndep.wr, by = "gid")
saveRDS(all.wr.sf, "WR/all.wr.sf.rds")The data was then cleaned and the centroid of each coastal drainage basin was computed.
all.wr.sf <- readRDS("WR/all.wr.sf.rds")
wr.sf.95 <- all.wr.sf %>% dplyr::filter(ndvi != 0) %>% filter(runoff > 0) %>% filter(is.na(bogs) == F) %>% filter(is.na(alt) == F) %>% filter(is.na(slope) == F)
wr.sf.95 <- wr.sf.95[!duplicated(wr.sf.95$gid),]
wr.sf.95$Runoff <- wr.sf.95$runoff *365*24*60*60 # from kg/m2/s to kg/m2/year or mm/y
wr.sf.95$logRunoff <- log(wr.sf.95$Runoff)
sfc_point_to_matrix = function(x) {
matrix(
unlist(x, use.names = FALSE),
nrow = length(x),
byrow = TRUE,
dimnames = list(1:length(x))
)
}
wr.centroid <- st_centroid(wr.sf.95, of.largest.polygon = T)
coord <- sfc_point_to_matrix(wr.centroid$geometry) %>% as.data.frame() %>% setNames(c("Longitude","Latitude"))
wr.sf.95 <- cbind(wr.sf.95, coord)
wr.sf.95$area <- st_area(wr.sf.95)
saveRDS(wr.sf.95, "WR/wr.sf.95.rds")noforest <- filter(wr.sf.95,forest == 0) #%>% dplyr::select(noforest,c("gid","area"))
#ggplot(noforest[which.max(noforest$area),])+geom_sf(aes(fill=as.factor(gid)))
weird.gid <- noforest$gid[which.max(noforest$area)]
weird.cat <- wr.sf.95 %>% filter(gid == weird.gid) %>% st_transform(CRS("EPSG:3035"))
corine <- raster::raster("Corine_Land_Cover/U2006_CLC2000_V2020_20u1.tif")
#weird.clc <- raster::extract(corine,weird.cat)
#saveRDS(weird.clc,"WR/clc/weird.clc.rds")
weird.clc <- readRDS("WR/clc/weird.clc.rds")
weird.tab <- sapply(weird.clc, function(x) tabulate(x,45))
weird.tabl.area <- weird.tab*prod(res(corine))
catchment.area <- colSums(weird.tabl.area)
weird.tab.prop <- sweep(weird.tabl.area,2,catchment.area, FUN = "/")
wr.sf.95$bogs[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[36]
wr.sf.95$forest[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[24]+weird.tab.prop[24]+weird.tab.prop[25]
wr.sf.95$arable[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[12]+weird.tab.prop[16]+weird.tab.prop[18]+weird.tab.prop[19]+weird.tab.prop[20]
wr.sf.95$glacier[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[34]
wr.sf.95$bare[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[31]+weird.tab.prop[32]+weird.tab.prop[33]
saveRDS(wr.sf.95,"WR/wr.sf.95.rds")wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
toreplace <- which(names(wr.sf.95) %in% c("ndvi","bogs","forest","arable","glacier","bare","temp","prec","alt","slope","tndep"))
names(wr.sf.95)[toreplace] <- c("NDVI","Bog","Forest","Arable","Glacier","Bare","Temp","Precip","Alt","Slope","TNdep")
wr.sf.95$area <- as.numeric(wr.sf.95$area)
saveRDS(wr.sf.95, "WR/wr.sf.95.rds")wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
m <- c("NDVI","Runoff","Precip","Temp","Bog","Forest","Arable","TNdep")
u <- c("NDVI on scale 0 to 1", "Mean Runoff, mm/y","Mean monthly precipitation,\n mm","Mean annual temperature,\n °C","Proportion of bogs","Proportion of Forest","Proportion of\narable land","Total Nitrogen Deposition (mg/m2)")
couleurs <- c("Greens","Blues","Blues","RdYlBu","YlOrBr","RdYlGn","RdYlGn","PuRd")
dir <- c(T,F,F,T,T,F,T,F)
for(i in m){
mysf <- wr.sf.95
map <- ggplot(mysf)+geom_sf(aes(fill = .data[[i]], col = .data[[i]]))+
scale_fill_continuous_divergingx(palette = "RdYlBu",
aesthetics = c("colour","fill"),
rev = dir[grep(i,m)],
name = paste(u[grep(i,m)],"\n(",i,")",sep=""))+
theme_minimal(base_size = 9)+
theme(legend.position = "bottom")
filename_i <- paste("WR/wr.map",i,"png",sep = ".")
ggsave(plot = map, filename = filename_i, dpi = "retina", width = 10, height = 12, units = "cm")
}
2 Model prediction 1995
We model the mean TOC concentration in all the water region in 1995 based on the simple SEM model of the previous section, using NDVI, bogs and runoff as predictors.
2.1 Neighbour matrix
The neighbor matrix of the coastal drainage basins polygons were computed using their centroids.
wr.centroid <- st_centroid(wr.sf.95, of_largest_polygon = T)
wr.neigh.set <- knearneigh(wr.centroid, k = 100)
wr.neigh.nb <- knn2nb(wr.neigh.set)
wr.kmat <- nb2listw(wr.neigh.nb)
saveRDS(wr.kmat,"WR/wr.kmat.rds")2.2 Model TOC
The TOC concentration in coastal drainage basins was modeled using the Spatial Error Linear Model with 5 predictors:
library(spatialreg)
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred <- predict(sem.fennoscandia.5,wr.sf.95,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
sem.pred.wr <- cbind(wr.sf.95,wr.pred$fit)
sem.pred.wr$TOC95 <- exp(sem.pred.wr$wr.pred.fit)
map.fit <- ggplot()+geom_sf(data = st_as_sf(sem.pred.wr), aes(fill = TOC95, col = TOC95))+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1,
name = "Modelled TOC 1995 (mg/L)")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
ggsave(plot = map.fit, filename = "WR/wr.toc.png",dpi = "retina", width = 6, height = 8, units = "cm")
map.log <- ggplot()+geom_sf(data = sem.pred.wr, aes(fill = wr.pred.fit, col = wr.pred.fit))+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1,
name = "Modelled log(TOC)")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
ggsave(plot = map.log, filename = "WR/wr.log.toc.png",dpi = "retina", width = 6, height = 8, units = "cm")
saveRDS(sem.pred.wr,"WR/sem.pred.wr.rds")knitr::include_graphics("WR/wr.log.toc.png")
Figure 6: Map of log(TOC) in coastal drainage basins for 1995
3 Forecast
3.1 Forecast of predictors
3.1.1 Forecast NDVI
3.1.1.1 Extract future temperature and precipitation
The forecast of the temperature and precipitation in the periods 2041-2060 and 2081-2100, for SSP 1-2.6 and SSP 3-7.0, were extracted based on the Future Worldclim dataset.
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.cnrm.temp <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp126_2041-2060.tif", band = 1)
ssp370.cnrm.temp <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp370_2041-2060.tif", band = 1)
ssp126.cnrm.temp.2050.crop <- raster::crop(ssp126.cnrm.temp,c(0,35,55,73))
png("worldclim_future/ssp126.temp.2050.png")
plot(ssp126.cnrm.temp.2050.crop)
dev.off()
ssp370.cnrm.temp.2050.crop <- raster::crop(ssp370.cnrm.temp,c(0,35,55,73))
png("worldclim_future/ssp370.temp.2050.png")
plot(ssp370.cnrm.temp.2050.crop)
dev.off()
ssp126.temp.wr <- raster::extract(ssp126.cnrm.temp,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm =T)
names(ssp126.temp.wr) <- c("gid", "temp126")
saveRDS(ssp126.temp.wr,"WR/ssp126.temp.wr.rds")
ssp370.temp.wr <- raster::extract(ssp370.cnrm.temp,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.temp.wr) <- c("gid", "temp370")
saveRDS(ssp370.temp.wr,"WR/ssp370.temp.wr.rds")wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.temp.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp126_2081-2100.tif", band = 1)
raster.ssp370.temp.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp370_2081-2100.tif", band = 1)
ssp126.cnrm.temp.2100.crop <- raster::crop(raster.ssp126.temp.2100,c(0,35,55,73))
png("worldclim_future/ssp126.temp.2100.png")
plot(ssp126.cnrm.temp.2100.crop)
dev.off()
ssp370.cnrm.temp.2100.crop <- raster::crop(raster.ssp370.temp.2100,c(0,35,55,73))
png("worldclim_future/ssp370.temp.2100.png")
plot(ssp370.cnrm.temp.2100.crop)
dev.off()
ssp126.temp.2100 <- raster::extract(raster.ssp126.temp.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm =T)
names(ssp126.temp.2100) <- c("gid", "temp126.2100")
saveRDS(ssp126.temp.2100,"WR/ssp126.temp.2100.rds")
ssp370.temp.2100 <- raster::extract(raster.ssp370.temp.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.temp.2100) <- c("gid", "temp370.2100")
saveRDS(ssp370.temp.2100,"WR/ssp370.temp.2100.rds")listgrob <- list("worldclim_future/ssp126.temp.2050.png", "worldclim_future/ssp370.temp.2050.png","worldclim_future/ssp126.temp.2100.png","worldclim_future/ssp370.temp.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("a) SSP126-2050","b) SSP370-2050","c) SSP126-2100","d) SSP370-2100"), label_size = 10)
Figure 7: Map of forecasted temperature (ºC) for 2050 and 2100, under SSP 1-2.6 and SSP 3-7.0
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.cnrm.prec <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp126_2041-2060.tif", band = 12)
ssp370.cnrm.prec <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp370_2041-2060.tif", band = 12)
ssp126.cnrm.prec.crop <- raster::crop(ssp126.cnrm.prec,c(0,35,55,73))
png("worldclim_future/ssp126.2050.prec.png")
plot(ssp126.cnrm.prec.crop)
dev.off()
ssp370.cnrm.prec.crop <- raster::crop(ssp370.cnrm.prec,c(0,35,55,73))
png("worldclim_future/ssp370.2050.prec.png")
plot(ssp370.cnrm.prec.crop)
dev.off()
ssp126.prec.wr <- raster::extract(ssp126.cnrm.prec,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp126.prec.wr) <- c("gid", "prec126")
saveRDS(ssp126.prec.wr,"WR/ssp126.prec.wr.rds")
ssp370.prec.wr <- raster::extract(ssp370.cnrm.prec,dplyr::select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.prec.wr) <- c("gid", "prec370")
saveRDS(ssp370.prec.wr,"WR/ssp370.prec.wr.rds")wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.prec.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp126_2081-2100.tif", band = 12)
raster.ssp370.prec.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp370_2081-2100.tif", band = 12)
ssp126.cnrm.prec.2100.crop <- raster::crop(raster.ssp126.prec.2100,c(0,35,55,73))
png("worldclim_future/ssp126.2100.prec.png")
plot(ssp126.cnrm.prec.2100.crop)
dev.off()
png("worldclim_future/ssp370.2100.prec.png")
ssp370.cnrm.prec.2100.crop <- raster::crop(raster.ssp370.prec.2100,c(0,35,55,73))
plot(ssp370.cnrm.prec.2100.crop)
ssp126.prec.2100 <- raster::extract(raster.ssp126.prec.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp126.prec.2100) <- c("gid", "prec126.2100")
saveRDS(ssp126.prec.2100,"WR/ssp126.prec.2100.rds")
ssp370.prec.2100 <- raster::extract(raster.ssp370.prec.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.prec.2100) <- c("gid", "prec370.2100")
saveRDS(ssp370.prec.2100,"WR/ssp370.prec.2100.rds")wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.temp.wr <- readRDS("WR/ssp126.temp.wr.rds")
ssp370.temp.wr <- readRDS("WR/ssp370.temp.wr.rds")
ssp126.prec.wr <- readRDS("WR/ssp126.prec.wr.rds")
ssp370.prec.wr <- readRDS("WR/ssp370.prec.wr.rds")
ssp126.temp.2100 <- readRDS("WR/ssp126.temp.2100.rds")
ssp370.temp.2100 <- readRDS("WR/ssp370.temp.2100.rds")
ssp126.prec.2100 <- readRDS("WR/ssp126.prec.2100.rds")
ssp370.prec.2100 <- readRDS("WR/ssp370.prec.2100.rds")
# SSP126
df.ssp126 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude")) %>% merge(ssp126.temp.wr, by = "gid") %>% merge(ssp126.prec.wr, by = "gid") %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","geometry"))
saveRDS(df.ssp126, "WR/df.ssp126.rds")
# SSP 370
df.ssp370 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude")) %>% merge(ssp370.temp.wr, by = "gid") %>% merge(ssp370.prec.wr, by = "gid") %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","geometry"))
saveRDS(df.ssp370, "WR/df.ssp370.rds")
# SSP126 2100
df.ssp126.2100 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude")) %>% merge(ssp126.temp.2100, by = "gid") %>% merge(ssp126.prec.2100, by = "gid") %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","geometry"))
saveRDS(df.ssp126.2100, "WR/df.ssp126.2100.rds")
# SSP 370 2100
df.ssp370.2100 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude")) %>% merge(ssp370.temp.2100, by = "gid") %>% merge(ssp370.prec.2100, by = "gid") %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","geometry"))
saveRDS(df.ssp370.2100, "WR/df.ssp370.2100.rds")
listgrob <- list("worldclim_future/ssp126.2050.prec.png", "worldclim_future/ssp370.2050.prec.png","worldclim_future/ssp126.2100.prec.png","worldclim_future/ssp370.2100.prec.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("a) SSP126-2050","b) SSP370-2050","c) SSP126-2100","d) SSP370-2100"), label_size = 10)
Figure 8: Map of forecasted precipitations (mm/s) for 2050 and 2100, under SSP 1-2.6 and SSP 3-7.0
3.1.1.2 NDVI model based on past temperature and precipitation values
We model NDVI based on temperature and precipitation, and use this model to predict the future NDVI in 2050. Choice of future data: ssp 126 as best scenario, ssp 370 as worse scenario (Hausfather 2019; CORDEX 2021). CRNM according to Raju and Kumar (2020).
Using a SELM model for NDVI, we notice that high residuals are observed in catchments where original NDVI is lower than 0,4. Myers-Smith et al. (2020) show that the relationship between NDVI and the Arctic biomass is not linear: under 0,4, an increase in NDVI reflects the vegetalization of bare ground, while over 0,4, an increase of NDVI reflects an increased productivity in established forests.
Several models were tested to model the NDVI: a linear model, a spatial error linear model, a quadratic model and a polynomial model.
library(betareg)
library(VGAM)
wr.kmat <- readRDS("WR/wr.kmat.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
par(mfrow = c(2,2))
#LM
summary(m1 <- lm(NDVI ~ Temp * log(Precip), data = wr.sf.95))
##
## Call:
## lm(formula = NDVI ~ Temp * log(Precip), data = wr.sf.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.65442 -0.05966 0.01396 0.08507 0.28204
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.933661 0.077110 25.08 <2e-16 ***
## Temp -0.223104 0.019800 -11.27 <2e-16 ***
## log(Precip) -0.200762 0.011751 -17.09 <2e-16 ***
## Temp:log(Precip) 0.037692 0.003035 12.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1202 on 1379 degrees of freedom
## Multiple R-squared: 0.3328, Adjusted R-squared: 0.3314
## F-statistic: 229.3 on 3 and 1379 DF, p-value: < 2.2e-16
summary(m1 <- step(m1))
## Start: AIC=-5855.14
## NDVI ~ Temp * log(Precip)
##
## Df Sum of Sq RSS AIC
## <none> 19.937 -5855.1
## - Temp:log(Precip) 1 2.2294 22.166 -5710.5
##
## Call:
## lm(formula = NDVI ~ Temp * log(Precip), data = wr.sf.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.65442 -0.05966 0.01396 0.08507 0.28204
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.933661 0.077110 25.08 <2e-16 ***
## Temp -0.223104 0.019800 -11.27 <2e-16 ***
## log(Precip) -0.200762 0.011751 -17.09 <2e-16 ***
## Temp:log(Precip) 0.037692 0.003035 12.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1202 on 1379 degrees of freedom
## Multiple R-squared: 0.3328, Adjusted R-squared: 0.3314
## F-statistic: 229.3 on 3 and 1379 DF, p-value: < 2.2e-16
plot(predict(m1), wr.sf.95$NDVI, pch = ".", main="a ) NDVI - LM", ylab="Observed NDVI", xlab="Predicted",abline(0,1))
text(x = 0.8, y = 0.4, label = paste("r2 =", round(summary(m1)$adj.r.squared,2), sep= ""))
# SELM
ndvi.model <- errorsarlm(formula = NDVI~Temp+log(Precip), data = wr.sf.95, listw = wr.kmat)
AIC <- 2*(3-1)-2*ndvi.model$LL
plot(ndvi.model$fitted.value, wr.sf.95$NDVI, pch = ".", main="b) NDVI - SELM", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.8, y = 0.5, label = paste("r2 =", round(cor(ndvi.model$fitted.values, wr.sf.95$NDVI)^2,2), sep= ""))
saveRDS(ndvi.model, "WR/selm.ndvi.rds")
# quadratic
summary(m2 <- lm(NDVI ~ (Temp + I(Temp^2)) * (log(Precip) + I(log(Precip)^2)), data = wr.sf.95))
##
## Call:
## lm(formula = NDVI ~ (Temp + I(Temp^2)) * (log(Precip) + I(log(Precip)^2)),
## data = wr.sf.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.58396 -0.03938 0.01550 0.05716 0.33305
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.962e+00 1.207e+00 -1.626 0.104212
## Temp 4.702e-01 6.420e-01 0.732 0.464027
## I(Temp^2) -2.305e-02 9.468e-02 -0.243 0.807662
## log(Precip) 9.996e-01 3.580e-01 2.793 0.005301 **
## I(log(Precip)^2) -9.158e-02 2.647e-02 -3.460 0.000558 ***
## Temp:log(Precip) -1.477e-01 1.903e-01 -0.776 0.437725
## Temp:I(log(Precip)^2) 1.287e-02 1.407e-02 0.914 0.360730
## I(Temp^2):log(Precip) 2.093e-03 2.800e-02 0.075 0.940438
## I(Temp^2):I(log(Precip)^2) 7.153e-05 2.066e-03 0.035 0.972391
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1027 on 1374 degrees of freedom
## Multiple R-squared: 0.5146, Adjusted R-squared: 0.5118
## F-statistic: 182.1 on 8 and 1374 DF, p-value: < 2.2e-16
summary(m2 <- step(m2))
## Start: AIC=-6285.14
## NDVI ~ (Temp + I(Temp^2)) * (log(Precip) + I(log(Precip)^2))
##
## Df Sum of Sq RSS AIC
## - I(Temp^2):I(log(Precip)^2) 1 0.0000126 14.504 -6287.1
## - I(Temp^2):log(Precip) 1 0.0000590 14.504 -6287.1
## - Temp:log(Precip) 1 0.0063610 14.510 -6286.5
## - Temp:I(log(Precip)^2) 1 0.0088239 14.513 -6286.3
## <none> 14.504 -6285.1
##
## Step: AIC=-6287.14
## NDVI ~ Temp + I(Temp^2) + log(Precip) + I(log(Precip)^2) + Temp:log(Precip) +
## Temp:I(log(Precip)^2) + I(Temp^2):log(Precip)
##
## Df Sum of Sq RSS AIC
## <none> 14.504 -6287.1
## - Temp:log(Precip) 1 0.041012 14.545 -6285.2
## - Temp:I(log(Precip)^2) 1 0.056474 14.560 -6283.8
## - I(Temp^2):log(Precip) 1 0.077693 14.582 -6281.7
##
## Call:
## lm(formula = NDVI ~ Temp + I(Temp^2) + log(Precip) + I(log(Precip)^2) +
## Temp:log(Precip) + Temp:I(log(Precip)^2) + I(Temp^2):log(Precip),
## data = wr.sf.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.58391 -0.03944 0.01542 0.05717 0.33340
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.977804 1.117506 -1.770 0.076976 .
## Temp 0.490437 0.265351 1.848 0.064780 .
## I(Temp^2) -0.026321 0.007294 -3.609 0.000319 ***
## log(Precip) 1.004431 0.330080 3.043 0.002387 **
## I(log(Precip)^2) -0.091938 0.024300 -3.783 0.000161 ***
## Temp:log(Precip) -0.153721 0.077960 -1.972 0.048833 *
## Temp:I(log(Precip)^2) 0.013311 0.005753 2.314 0.020824 *
## I(Temp^2):log(Precip) 0.003061 0.001128 2.714 0.006732 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1027 on 1375 degrees of freedom
## Multiple R-squared: 0.5146, Adjusted R-squared: 0.5122
## F-statistic: 208.3 on 7 and 1375 DF, p-value: < 2.2e-16
plot(predict(m2), wr.sf.95$NDVI, pch = ".", main="c) NDVI - quadratic model", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.8, y = 0.4, label = paste("r2 =", round(summary(m2)$adj.r.squared,2), sep= ""))
# poly
summary(m3 <- lm(NDVI ~ poly(cbind(Temp, log(Precip)), degree = 5), data = wr.sf.95))
##
## Call:
## lm(formula = NDVI ~ poly(cbind(Temp, log(Precip)), degree = 5),
## data = wr.sf.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.57595 -0.04017 0.01037 0.05052 0.33179
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.704752 0.007771 90.693
## poly(cbind(Temp, log(Precip)), degree = 5)1.0 2.218061 0.491375 4.514
## poly(cbind(Temp, log(Precip)), degree = 5)2.0 -1.199191 0.584234 -2.053
## poly(cbind(Temp, log(Precip)), degree = 5)3.0 -1.234325 0.380678 -3.242
## poly(cbind(Temp, log(Precip)), degree = 5)4.0 0.172241 0.233386 0.738
## poly(cbind(Temp, log(Precip)), degree = 5)5.0 0.074829 0.116072 0.645
## poly(cbind(Temp, log(Precip)), degree = 5)0.1 -1.193904 0.442951 -2.695
## poly(cbind(Temp, log(Precip)), degree = 5)1.1 -16.692855 30.793091 -0.542
## poly(cbind(Temp, log(Precip)), degree = 5)2.1 116.680490 34.912996 3.342
## poly(cbind(Temp, log(Precip)), degree = 5)3.1 -38.158776 21.683668 -1.760
## poly(cbind(Temp, log(Precip)), degree = 5)4.1 -1.452680 11.647690 -0.125
## poly(cbind(Temp, log(Precip)), degree = 5)0.2 1.181823 0.294145 4.018
## poly(cbind(Temp, log(Precip)), degree = 5)1.2 -91.265780 17.942551 -5.087
## poly(cbind(Temp, log(Precip)), degree = 5)2.2 83.064904 17.238786 4.818
## poly(cbind(Temp, log(Precip)), degree = 5)3.2 -56.596202 11.011082 -5.140
## poly(cbind(Temp, log(Precip)), degree = 5)0.3 0.708308 0.183026 3.870
## poly(cbind(Temp, log(Precip)), degree = 5)1.3 -58.065996 9.442858 -6.149
## poly(cbind(Temp, log(Precip)), degree = 5)2.3 34.180342 8.528169 4.008
## poly(cbind(Temp, log(Precip)), degree = 5)0.4 0.142940 0.125775 1.136
## poly(cbind(Temp, log(Precip)), degree = 5)1.4 2.984816 6.343738 0.471
## poly(cbind(Temp, log(Precip)), degree = 5)0.5 -0.081876 0.108766 -0.753
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## poly(cbind(Temp, log(Precip)), degree = 5)1.0 6.91e-06 ***
## poly(cbind(Temp, log(Precip)), degree = 5)2.0 0.040303 *
## poly(cbind(Temp, log(Precip)), degree = 5)3.0 0.001214 **
## poly(cbind(Temp, log(Precip)), degree = 5)4.0 0.460637
## poly(cbind(Temp, log(Precip)), degree = 5)5.0 0.519247
## poly(cbind(Temp, log(Precip)), degree = 5)0.1 0.007118 **
## poly(cbind(Temp, log(Precip)), degree = 5)1.1 0.587840
## poly(cbind(Temp, log(Precip)), degree = 5)2.1 0.000854 ***
## poly(cbind(Temp, log(Precip)), degree = 5)3.1 0.078667 .
## poly(cbind(Temp, log(Precip)), degree = 5)4.1 0.900765
## poly(cbind(Temp, log(Precip)), degree = 5)0.2 6.19e-05 ***
## poly(cbind(Temp, log(Precip)), degree = 5)1.2 4.15e-07 ***
## poly(cbind(Temp, log(Precip)), degree = 5)2.2 1.61e-06 ***
## poly(cbind(Temp, log(Precip)), degree = 5)3.2 3.15e-07 ***
## poly(cbind(Temp, log(Precip)), degree = 5)0.3 0.000114 ***
## poly(cbind(Temp, log(Precip)), degree = 5)1.3 1.02e-09 ***
## poly(cbind(Temp, log(Precip)), degree = 5)2.3 6.46e-05 ***
## poly(cbind(Temp, log(Precip)), degree = 5)0.4 0.255959
## poly(cbind(Temp, log(Precip)), degree = 5)1.4 0.638063
## poly(cbind(Temp, log(Precip)), degree = 5)0.5 0.451717
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09847 on 1362 degrees of freedom
## Multiple R-squared: 0.5581, Adjusted R-squared: 0.5516
## F-statistic: 86 on 20 and 1362 DF, p-value: < 2.2e-16
summary(m3 <- step(m3))
## Start: AIC=-6390.86
## NDVI ~ poly(cbind(Temp, log(Precip)), degree = 5)
##
## Df Sum of Sq RSS AIC
## <none> 13.206 -6390.9
## - poly(cbind(Temp, log(Precip)), degree = 5) 20 16.677 29.882 -5301.5
##
## Call:
## lm(formula = NDVI ~ poly(cbind(Temp, log(Precip)), degree = 5),
## data = wr.sf.95)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.57595 -0.04017 0.01037 0.05052 0.33179
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.704752 0.007771 90.693
## poly(cbind(Temp, log(Precip)), degree = 5)1.0 2.218061 0.491375 4.514
## poly(cbind(Temp, log(Precip)), degree = 5)2.0 -1.199191 0.584234 -2.053
## poly(cbind(Temp, log(Precip)), degree = 5)3.0 -1.234325 0.380678 -3.242
## poly(cbind(Temp, log(Precip)), degree = 5)4.0 0.172241 0.233386 0.738
## poly(cbind(Temp, log(Precip)), degree = 5)5.0 0.074829 0.116072 0.645
## poly(cbind(Temp, log(Precip)), degree = 5)0.1 -1.193904 0.442951 -2.695
## poly(cbind(Temp, log(Precip)), degree = 5)1.1 -16.692855 30.793091 -0.542
## poly(cbind(Temp, log(Precip)), degree = 5)2.1 116.680490 34.912996 3.342
## poly(cbind(Temp, log(Precip)), degree = 5)3.1 -38.158776 21.683668 -1.760
## poly(cbind(Temp, log(Precip)), degree = 5)4.1 -1.452680 11.647690 -0.125
## poly(cbind(Temp, log(Precip)), degree = 5)0.2 1.181823 0.294145 4.018
## poly(cbind(Temp, log(Precip)), degree = 5)1.2 -91.265780 17.942551 -5.087
## poly(cbind(Temp, log(Precip)), degree = 5)2.2 83.064904 17.238786 4.818
## poly(cbind(Temp, log(Precip)), degree = 5)3.2 -56.596202 11.011082 -5.140
## poly(cbind(Temp, log(Precip)), degree = 5)0.3 0.708308 0.183026 3.870
## poly(cbind(Temp, log(Precip)), degree = 5)1.3 -58.065996 9.442858 -6.149
## poly(cbind(Temp, log(Precip)), degree = 5)2.3 34.180342 8.528169 4.008
## poly(cbind(Temp, log(Precip)), degree = 5)0.4 0.142940 0.125775 1.136
## poly(cbind(Temp, log(Precip)), degree = 5)1.4 2.984816 6.343738 0.471
## poly(cbind(Temp, log(Precip)), degree = 5)0.5 -0.081876 0.108766 -0.753
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## poly(cbind(Temp, log(Precip)), degree = 5)1.0 6.91e-06 ***
## poly(cbind(Temp, log(Precip)), degree = 5)2.0 0.040303 *
## poly(cbind(Temp, log(Precip)), degree = 5)3.0 0.001214 **
## poly(cbind(Temp, log(Precip)), degree = 5)4.0 0.460637
## poly(cbind(Temp, log(Precip)), degree = 5)5.0 0.519247
## poly(cbind(Temp, log(Precip)), degree = 5)0.1 0.007118 **
## poly(cbind(Temp, log(Precip)), degree = 5)1.1 0.587840
## poly(cbind(Temp, log(Precip)), degree = 5)2.1 0.000854 ***
## poly(cbind(Temp, log(Precip)), degree = 5)3.1 0.078667 .
## poly(cbind(Temp, log(Precip)), degree = 5)4.1 0.900765
## poly(cbind(Temp, log(Precip)), degree = 5)0.2 6.19e-05 ***
## poly(cbind(Temp, log(Precip)), degree = 5)1.2 4.15e-07 ***
## poly(cbind(Temp, log(Precip)), degree = 5)2.2 1.61e-06 ***
## poly(cbind(Temp, log(Precip)), degree = 5)3.2 3.15e-07 ***
## poly(cbind(Temp, log(Precip)), degree = 5)0.3 0.000114 ***
## poly(cbind(Temp, log(Precip)), degree = 5)1.3 1.02e-09 ***
## poly(cbind(Temp, log(Precip)), degree = 5)2.3 6.46e-05 ***
## poly(cbind(Temp, log(Precip)), degree = 5)0.4 0.255959
## poly(cbind(Temp, log(Precip)), degree = 5)1.4 0.638063
## poly(cbind(Temp, log(Precip)), degree = 5)0.5 0.451717
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09847 on 1362 degrees of freedom
## Multiple R-squared: 0.5581, Adjusted R-squared: 0.5516
## F-statistic: 86 on 20 and 1362 DF, p-value: < 2.2e-16
plot(predict(m3), wr.sf.95$NDVI, pch = ".", main="d) NDVI - polynomial ", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.8, y = 0.4, label = paste("r2 =", round(summary(m3)$adj.r.squared,2), sep= ""))
saveRDS(m3,"WR/poly.ndvi.rds")
Figure 9: Test of 4 NDVI models, a) linear model, b) spatial error linear model, c) quadratic model, d) polynomial model
The SELM and the polynomial model had the best results. We tested them on the Fennoscandia dataset.
k.neigh.w <- readRDS("k.neigh.w.rds")
fennoscandia <- readRDS("fennoscandia.rds")
fennoscandia.33 <- readRDS("fennoscandia.33.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
nor <- st_read("Country_shapefile/norway.shp")
swe <- st_read("Country_shapefile/sweden.shp")
fin <- st_read("Country_shapefile/finland.shp")
selm.ndvi <- readRDS("WR/selm.ndvi.rds")
selm.ndvi.fen <- predict(selm.ndvi,fennoscandia.33,k.neigh.w, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
fennoscandia.selm.ndvi.model <- fennoscandia.33 %>% st_as_sf() %>% mutate(ndvi.fit = selm.ndvi.fen$fit, ndvi.residuals = fennoscandia.33$NDVI - selm.ndvi.fen$fit)
fennoscandia.selm.ndvi.model$area <- st_area(fennoscandia.selm.ndvi.model) %>% as.numeric()
ggplot(fennoscandia.selm.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(fennoscandia.selm.ndvi.model$NDVI,fennoscandia.selm.ndvi.model$ndvi.fit)^2,2)), x = 0.75, y = 0.3, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/fen.ndvi.selm.results.png", dpi = "retina", width = 6, height = 8, units = "cm")
fen.residuals.selm <- ggplot()+
geom_sf(data = fennoscandia.selm.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "SELM residuals")+
theme_minimal(base_size = 6)+theme(legend.position = "bottom")
ggsave(plot = fen.residuals.selm, filename = "WR/fen.selm.ndvi.residuals.png",dpi = "retina", width = 6, height = 8, units = "cm")
# poly.ndvi
poly.ndvi <- readRDS("WR/poly.ndvi.rds")
poly.ndvi.fen <- predict(poly.ndvi,fennoscandia.33)
fennoscandia.poly.ndvi.model <- fennoscandia.33 %>% mutate(ndvi.fit = poly.ndvi.fen, ndvi.residuals = fennoscandia.33$NDVI - poly.ndvi.fen)
fennoscandia.poly.ndvi.model$area <- st_area(fennoscandia.selm.ndvi.model) %>% as.numeric()
resi <- filter(select(fennoscandia.poly.ndvi.model, c(ebint, nation, area, TOC, NDVI,ndvi.residuals)), ndvi.residuals < -0.2)
ggplot(resi)+geom_sf(aes(col = ebint))
test.compare <- filter(fennoscandia.poly.ndvi.model, between(ebint, 16600000, 16620000)) %>% select(c("ebint","nation","TOC","NDVI","ndvi.residuals","Runoff","Temp","Precip","Slope","Bog","Forest","Arable","TNdep"))
ggplot(fennoscandia.poly.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(fennoscandia.poly.ndvi.model$NDVI,fennoscandia.poly.ndvi.model$ndvi.fit),2)), x = 0.75, y = 0.3, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/fen.ndvi.poly.results.png", dpi = "retina", width = 6, height = 8, units = "cm")
fen.residuals <- ggplot()+
geom_sf(data = fennoscandia.poly.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "Polynomial residuals")+
theme_minimal(base_size = 6)+theme(legend.position = "bottom")
ggsave(plot = fen.residuals, filename = "WR/fen.poly.ndvi.residuals.png",dpi = "retina", width = 6, height = 8, units = "cm")
listgrob <- list("WR/fen.ndvi.selm.results.png","WR/fen.ndvi.poly.results.png","WR/fen.selm.ndvi.residuals.png","WR/fen.poly.ndvi.residuals.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol = 2, nrow = 2, rel_heights = c(1,2), labels = c("a) Observed vs predicted NDVI for SELM","b) Observed vs predicted NDVI for polynomial model","c) Map of residuals for SELM","d) Map of residuals for Polynomial model"), hjust = 3)
Figure 10: Results of NDVI models (SELM and polynomial): observed vs predicted for a) SELM and b) polynoimal model, and map of residuals (observed - predicted) for c) SELM and d) Polynomial
The polynomial model preforms better on the test set (Fennoscandia catchments). This model will be used to forecast NDVI in 2050 and in 2100.
3.1.1.3 Forecast NDVI in 2050 and 2100
df.ssp126 <- readRDS("WR/df.ssp126.rds")
df.ssp370 <- readRDS("WR/df.ssp370.rds")
df.ssp126.2100 <- readRDS("WR/df.ssp126.2100.rds")
df.ssp370.2100 <- readRDS("WR/df.ssp370.2100.rds")
ndvi.model <- readRDS("WR/ndvi.model.rds")
wr.kmat <- readRDS("WR/wr.kmat.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
#NDVI 126 2050
ndvi.ssp126.fit <- predict(object = ndvi.model, newdata = df.ssp126, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
ndvi.ssp126 <- cbind(df.ssp126, ndvi.ssp126.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
saveRDS(ndvi.ssp126,"WR/ndvi.ssp126.rds")
# NDVI 370 2050
ndvi.ssp370.fit <- predict(object = ndvi.model, newdata = df.ssp370, listw = wr.kmat, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
ndvi.ssp370 <- cbind(df.ssp370, ndvi.ssp370.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
saveRDS(ndvi.ssp370,"WR/ndvi.ssp370.rds")
#NDVI 126 2010
ndvi.ssp126.2100.fit <- predict(object = ndvi.model, newdata = df.ssp126.2100, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
ndvi.ssp126.2100 <- cbind(df.ssp126.2100, ndvi.ssp126.2100.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
saveRDS(ndvi.ssp126.2100,"WR/ndvi.ssp126.2100.rds")
#NDVI 370 2100
ndvi.ssp370.2100.fit <- predict(object = ndvi.model, newdata = df.ssp370.2100, listw = wr.kmat, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% setNames(c("ndvi","trend","signal"))
ndvi.ssp370.2100 <- cbind(df.ssp370.2100, ndvi.ssp370.2100.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
saveRDS(ndvi.ssp370.2100,"WR/ndvi.ssp370.2100.rds")
Figure 10: Results of NDVI models (SELM and polynomial): observed vs predicted for a) SELM and b) polynoimal model, and map of residuals (observed - predicted) for c) SELM and d) Polynomial
3.1.2 Forecast runoff 2050 and 2100
In the historical dataset in Worldclim, the runoff is given as a yearly average in kg/m2/s, or mm/s.
In the CNRM-CM6 forecast we get the monthly forecast for runoff (in mm/s) and convert it in mm/y.
listgrob <- list("CORDEX_future/runoff.ssp126.2050.png", "CORDEX_future/runoff.ssp370.2050.png","CORDEX_future/runoff.ssp126.2100.png","CORDEX_future/runoff.ssp370.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
ggarrange(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("SSP126-2050","SSP370-2050","SSP126-2100","SSP370-2100"))
wr.sf <- readRDS("WR/wr.sf.rds")
ssp126.runoff.wr <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp126_r1i1p1f2_gr_201501-210012.nc", bands = c(313:552))
ssp126.runoff.wr.mean <- ssp126.runoff.wr %>% mean()
ssp126.runoff.wr.mean.crop <- raster::crop(ssp126.runoff.wr.mean,c(0,35,55,73))
plot(ssp126.runoff.wr.mean.crop)
ssp126.runoff <- raster::extract(ssp126.runoff.wr.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp126.runoff) <- c("gid","runoff")
ssp126.runoff$Runoff <- ssp126.runoff$runoff *365*24*60*60
saveRDS(ssp126.runoff,"WR/ssp126.runoff.2050.rds")
ssp370.runoff.wr <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp370_r1i1p1f2_gr_201501-210012.nc", bands = c(313:552))
ssp370.runoff.wr.mean <- ssp370.runoff.wr %>% mean()
ssp370.runoff.wr.mean.crop <- raster::crop(ssp370.runoff.wr.mean,c(0,35,55,73))
plot(ssp370.runoff.wr.mean.crop)
ssp370.runoff <- raster::extract(ssp370.runoff.wr.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp370.runoff) <- c("gid","runoff")
ssp370.runoff$Runoff <- ssp370.runoff$runoff *365*24*60*60
saveRDS(ssp370.runoff,"WR/ssp370.runoff.2050.rds")wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.runoff.2100 <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp126_r1i1p1f2_gr_201501-210012.nc", bands = c(805:1032))
raster.ssp126.runoff.2100.mean <- raster.ssp126.runoff.2100 %>% mean()
raster.ssp126.runoff.2100.mean.crop <- raster::crop(raster.ssp126.runoff.2100.mean,c(0,35,55,73))
plot(raster.ssp126.runoff.2100.mean.crop)
ssp126.runoff.2100 <- raster::extract(raster.ssp126.runoff.2100.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp126.runoff.2100) <- c("gid","runoff")
ssp126.runoff.2100$Runoff <- ssp126.runoff.2100$runoff *365*24*60*60
saveRDS(ssp126.runoff.2100,"WR/ssp126.runoff.2100.rds")
raster.ssp370.runoff.2100 <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp370_r1i1p1f2_gr_201501-210012.nc", bands = c(805:1032))
raster.ssp370.runoff.2100.mean <- raster.ssp370.runoff.2100 %>% mean()
raster.ssp370.runoff.2100.mean.crop <- raster::crop(raster.ssp370.runoff.2100.mean,c(0,35,55,73))
plot(raster.ssp370.runoff.2100.mean.crop)
ssp370.runoff.2100 <- raster::extract(raster.ssp370.runoff.2100.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp370.runoff.2100) <- c("gid","runoff")
ssp370.runoff.2100$Runoff <- ssp370.runoff.2100$runoff *365*24*60*60
saveRDS(ssp370.runoff.2100,"WR/ssp370.runoff.2100.rds")3.1.3 Future TNdep
Future TNdep will largely depend on the public policies implemented under the scenarios SSP 1-2.6 and 3-7.0.
TNdep data are from 2000.
Emission reduction between 2005 and 2020 (European Environment Agency 2021):
- NH3: -8%
- NOx -42%
Predictions for future emissions for NOx (summary for policymakers AR6 (Grubb et al., n.d.)):
- SSP1-2.6: from 11 MtN2O/y to 9 MtN2O (18% decrease) in 2050, stable in 2100 (also 20% compared to 2015)
- SSP3-7.0: from 11 MtN2o/y to 15 MtN2O (137% increase) in 2050, and 21 MtN2O (190% increase) in 2100
Assumption for NH3:
- SSP1-2.6: Decrease according to targets (-20%)
- SSP3-7.0: No decrease
ndep.wr <- readRDS("WR/ndep.wr.rds")
ndep.sum <- ndep.wr@data %>% apply(2, sum)
(ndep.sum["woxn"]+ndep.sum["doxn"])/ndep.sum["tndep"]
(ndep.sum["wrdn"]+ndep.sum["drdn"])/ndep.sum["tndep"]
# SSP126 2050
ssp126.ndep.2050 <- data.frame(gid = ndep.wr$gid)
ssp126.ndep.2050$oxn <- ndep.wr$woxn*0.55*0.8 + ndep.wr$doxn*0.55*0.8 # -45% for decrease between 2000 and 2020, -20% between 2020 and 2050
ssp126.ndep.2050$rdn <- ndep.wr$wrdn*0.9*0.8 + ndep.wr$drdn*0.9*0.8 # -10% for decrease between 2000 and 2020, -20% between 2020 and 2050
ssp126.ndep.2050$TNdep <- ssp126.ndep.2050$oxn + ssp126.ndep.2050$rdn
saveRDS(ssp126.ndep.2050,"WR/ssp126.ndep.2050.rds")
# SSP126 2100
ssp126.ndep.2100 <- data.frame(gid = ndep.wr$gid)
ssp126.ndep.2100$oxn <- ndep.wr$woxn*0.55*0.8+ndep.wr$doxn*0.55*0.8 # -45% for decrease between 2000 and 2020, -20% between 2020 and 2100
ssp126.ndep.2100$rdn <- ndep.wr$wrdn*0.9*0.8+ndep.wr$drdn*0.9*0.8 # -10% for decrease between 2000 and 2020, -20% between 2020 and 2100
ssp126.ndep.2100$TNdep <- ssp126.ndep.2100$oxn + ssp126.ndep.2100$rdn
saveRDS(ssp126.ndep.2100,"WR/ssp126.ndep.2100.rds")
# SSP370 2050
ssp370.ndep.2050 <- data.frame(gid = ndep.wr$gid)
ssp370.ndep.2050$oxn <- ndep.wr$woxn*0.55*1.37+ndep.wr$doxn*0.55*1.37 # -45% for decrease between 2000 and 2020, 137% increase between 2020 and 2100
ssp370.ndep.2050$rdn <- ndep.wr$wrdn*0.9+ndep.wr$drdn*0.9 # -10% for decrease between 2000 and 2020, no change between 2020 and 2100
ssp370.ndep.2050$TNdep <- ssp370.ndep.2050$oxn + ssp370.ndep.2050$rdn
saveRDS(ssp370.ndep.2050,"WR/ssp370.ndep.2050.rds")
# SSP370 2100
ssp370.ndep.2100 <- data.frame(gid = ndep.wr$gid)
ssp370.ndep.2100$oxn <- ndep.wr$woxn*0.55*1.9+ndep.wr$doxn*0.55*1.9 # -45% for decrease between 2000 and 2020, 190% increase between 2020 and 2100
ssp370.ndep.2100$rdn <- ndep.wr$wrdn*0.9+ndep.wr$drdn*0.9 # -10% for decrease between 2000 and 2020, no change between 2020 and 2100
ssp370.ndep.2100$TNdep <- ssp370.ndep.2100$oxn + ssp370.ndep.2100$rdn
saveRDS(ssp370.ndep.2100,"WR/ssp370.ndep.2100.rds")3.1.4 Gather datasets and plot
The datasets were gathered into a single dataframe.
ndvi.ssp126 <- readRDS("WR/ndvi.ssp126.rds")
ndvi.ssp370 <- readRDS("WR/ndvi.ssp370.rds")
ssp126.runoff.2050 <- readRDS("WR/ssp126.runoff.2050.rds")
ssp370.runoff.2050 <- readRDS("WR/ssp370.runoff.2050.rds")
ssp126.ndep.2050 <- readRDS("WR/ssp126.ndep.2050.rds")
ssp370.ndep.2050 <- readRDS("WR/ssp370.ndep.2050.rds")
ndvi.ssp126.2100 <- readRDS("WR/ndvi.ssp126.2100.rds")
ndvi.ssp370.2100 <- readRDS("WR/ndvi.ssp370.2100.rds")
ssp126.runoff.2100 <- readRDS("WR/ssp126.runoff.2100.rds")
ssp370.runoff.2100 <- readRDS("WR/ssp370.runoff.2100.rds")
ssp126.ndep.2100 <- readRDS("WR/ssp126.ndep.2100.rds")
ssp370.ndep.2100 <- readRDS("WR/ssp370.ndep.2100.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
# 2050 -----
wr.sf.126.2050 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(st_drop_geometry(ndvi.ssp126), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp126.runoff.2050)), by = "gid") %>% merge(ssp126.ndep.2050)
wr.sf.370.2050 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(st_drop_geometry(ndvi.ssp370), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp370.runoff.2050)), by = "gid") %>% merge(ssp370.ndep.2050)
na.lines.126 <- which(complete.cases(st_drop_geometry(wr.sf.126.2050))==F)
na.lines.370 <- which(complete.cases(st_drop_geometry(wr.sf.370.2050))==F)
wr.sf.126.2050 <- wr.sf.126.2050[-na.lines.126,]
wr.sf.370.2050 <- wr.sf.370.2050[-na.lines.370,]
wr.sf.126.2050$logRunoff <- wr.sf.126.2050$Runoff %>% log()
wr.sf.370.2050$logRunoff <- wr.sf.370.2050$Runoff %>% log()
saveRDS(wr.sf.126.2050, "WR/wr.sf.126.2050.rds")
saveRDS(wr.sf.370.2050, "WR/wr.sf.370.2050.rds")
# 2100 ----
wr.sf.126.2100 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(st_drop_geometry(ndvi.ssp126.2100), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp126.runoff.2100)), by = "gid") %>% merge(ssp126.ndep.2100)
wr.sf.370.2100 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(st_drop_geometry(ndvi.ssp370.2100), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp370.runoff.2100)), by = "gid") %>% merge(ssp370.ndep.2100)
na.lines.126 <- which(complete.cases(st_drop_geometry(wr.sf.126.2100))==F)
na.lines.370 <- which(complete.cases(st_drop_geometry(wr.sf.370.2100))==F)
wr.sf.126.2100 <- wr.sf.126.2100[-na.lines.126,]
wr.sf.370.2100 <- wr.sf.370.2100[-na.lines.370,]
wr.sf.126.2100$logRunoff <- wr.sf.126.2100$Runoff %>% log()
wr.sf.370.2100$logRunoff <- wr.sf.370.2100$Runoff %>% log()
saveRDS(wr.sf.126.2100, "WR/wr.sf.126.2100.rds")
saveRDS(wr.sf.370.2100, "WR/wr.sf.370.2100.rds")
Then we plot the difference between the forecasted variables in 2050 and 2100, compared to 1995.
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
diff.var <- function(sf, vars1, vars2){
sf2 <- sf
for (i in 1:length(vars2)){
newname <- paste("diff", vars2[i], sep = "")
if(vars2[i] %in% c("Temp","NDVI","TNdep")){
x <- sf2[,vars2[i]] %>% st_drop_geometry()
y <- sf2[,vars1[i]] %>% st_drop_geometry()
z <- x-y
names(z) <- newname
sf2 <- cbind(sf2, z)
} else if(vars2[i] %in% c("Precip", "Runoff")){
x <- sf2[,vars2[i]] %>% st_drop_geometry()
y <- sf2[,vars1[i]] %>% st_drop_geometry()
z <- x - y
names(z) <- newname
sf2 <- cbind(sf2, z)
}
}
return(sf2)
}
wr.sf.126.2050 <- diff.var(wr.sf.126.2050, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.370.2050 <- diff.var(wr.sf.370.2050, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.126.2100 <- diff.var(wr.sf.126.2100, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.370.2100 <- diff.var(wr.sf.370.2100, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
list.sf <- list(wr.sf.126.2050, wr.sf.126.2100, wr.sf.370.2050, wr.sf.370.2100)
names(list.sf) <- c("SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
all.diff <- data.frame(gid = c(wr.sf.126.2050$gid, wr.sf.370.2050$gid, wr.sf.126.2100$gid, wr.sf.370.2100$gid))
for(x in c("diffTemp","diffPrecip","diffNDVI","diffRunoff","diffTNdep")){
all.vec <- c()
for(i in 1:length(list.sf)){
vec <- select(st_drop_geometry(list.sf[[i]]),all_of(x))
all.vec <<- rbind(all.vec,vec)
}
all.diff <- cbind(all.diff, all.vec)
}
max.lim <- apply(all.diff, 2, max, na.rm = T)
min.lim <- apply(all.diff, 2, min, na.rm = T)
# x1 <- filter(wr.sf.126.2050, NDVI95 > 0.4)
# x2 <- filter(wr.sf.370.2050, NDVI95 > 0.4)
# x3 <- filter(wr.sf.126.2100, NDVI95 > 0.4)
# x4 <- filter(wr.sf.370.2100, NDVI95 > 0.4)
# all.ndvi <- c(x1$NDVI/x1$NDVI95*100-100,x2$NDVI/x2$NDVI95*100-100, x3$NDVI/x3$NDVI95*100-100, x4$NDVI/x4$NDVI95*100-100)
# palettes <- data.frame(diffTemp = brewer.pal(name = "RdYlBu", n = 9)[c(1,5,9)],
# diffPrecip = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)],
# diffNDVI = brewer.pal(name = "RdYlGn", n = 9)[c(9,5,1)],
# diffRunoff = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)],
# diffTNdep = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)])
legends <- c("Difference Temperature (°C)", "Difference Precipitation (mm/y)", "Difference NDVI", "Difference Runoff (mm/y)","DIfference Nitrogen Deposition (ug/m2)")
diffs <- c("diffTemp","diffPrecip","diffNDVI","diffRunoff","diffTNdep")
rev.pal <- c(T,F,T,F,T)
for(i in c(1:length(list.sf))){
for(x in diffs){
g <- ggplot(list.sf[[i]])+geom_sf(aes(fill = .data[[x]], col = .data[[x]])) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = rev.pal[grep(x, diffs)] , aesthetics = c("fill","col"), name = paste(legends[grep(x, diffs)],"\n",names(list.sf)[i], sep = ""), limits = c(min.lim[[x]], max.lim[[x]])) +
# scale_fill_gradient2(low=palettes[3, grep(x,names(palettes))],mid=palettes[2, grep(x,names(palettes))],high = palettes[1, grep(x,names(palettes))],
# midpoint = 0, aesthetics = c("fill","col"),
# name = paste(legends[grep(x, diffs)],"\n",names(list.sf)[i], sep = ""), limits = c(min.lim[[x]], max.lim[[x]])) +
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = g, filename = paste("WR/", names(list.sf)[i], ".", x, ".png",sep = ""),dpi = "retina", width = 12, height = 16, units = "cm")
}
}
white <- ggplot()+theme_void()
ggsave(plot = white, filename = "white.png",dpi = "retina", width = 12, height = 16, units = "cm")listgrob <- list("WR/SSP126.2050.diffTemp.png", "WR/SSP126.2100.diffTemp.png","WR/SSP370.2050.diffTemp.png","WR/SSP370.2100.diffTemp.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
Figure 7: Map of forecasted temperature (ºC) for 2050 and 2100, under SSP 1-2.6 and SSP 3-7.0
listgrob <- list("WR/SSP126.2050.diffPrecip.png","WR/SSP126.2100.diffPrecip.png","WR/SSP370.2050.diffPrecip.png","WR/SSP370.2100.diffPrecip.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
Figure 11: Forecast of mean precipitation intensities for the period 2041-2060 (i.e., 2050) and for the period 2081-2100 (i.e., 2100), and change in precipitation compared to the average for 1995 to 2050, under SSP 1-2.6 (2ºC target) and SSP 3-7.0 (realistic worst-case) scenarios.
listgrob <- list("WR/SSP126.2050.diffNDVI.png","WR/SSP126.2100.diffNDVI.png","WR/SSP370.2050.diffNDVI.png","WR/SSP370.2100.diffNDVI.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
Figure 12: Forecast of mean NDVI for the period 2041-2060 (i.e., 2050) and for the period 2081-2100 (i.e., 2100), and change in NDVI compared to the average for 1995 to 2050, under SSP 1-2.6 (2ºC target) and SSP 3-7.0 (realistic worst-case) scenarios.
listgrob <- list("WR/SSP126.2050.diffRunoff.png","WR/SSP126.2100.diffRunoff.png","WR/SSP370.2050.diffRunoff.png","WR/SSP370.2100.diffRunoff.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
Figure 13: Forecast of mean runoff intensities for the period 2041-2060 (i.e., 2050) and for the period 2081-2100 (i.e., 2100), and change in runoff compared to the average for 1995 to 2050, under SSP 1-2.6 (2ºC target) and SSP 3-7.0 (realistic worst-case) scenarios.
listgrob <- list("WR/SSP126.2050.diffTNdep.png","WR/SSP126.2100.diffTNdep.png","WR/SSP370.2050.diffTNdep.png","WR/SSP370.2100.diffTNdep.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
Figure 14: Forecast of mean total nitrogen deposition for the period 2041-2060 (i.e., 2050) and for the period 2081-2100 (i.e., 2100), and change in TNdep compared to the average for 1995 to 2050, under SSP 1-2.6 (2ºC target) and SSP 3-7.0 (realistic worst-case) scenarios.
3.2 Forecast TOC for 2050 and 2100
The average concentration of total organic carbon in the coastal drainage basins was predicted based on the SELM with 5 predictors, fitted on the Fennoscandian dataset. The difference between the predicted log(TOC) for 2050 and 2100 under the SSP 1-2.6 and 3-7.0 scenarios, compared to the log(TOC) predicted for 1995, was computed and plotted.
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp126 <- predict(sem.fennoscandia.5,wr.sf.126.2050,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp126 <- cbind(wr.sf.126.2050,wr.pred.ssp126)
wr.ssp126$TOC <- exp(wr.ssp126$fit)
wr.ssp126 <- merge(wr.ssp126, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp126,"WR/wr.ssp126.2050.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp370 <- predict(sem.fennoscandia.5,wr.sf.370.2050,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp370 <- cbind(wr.sf.370.2050,wr.pred.ssp370)
wr.ssp370$TOC <- exp(wr.ssp370$fit)
wr.ssp370 <- merge(wr.ssp370, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp370,"WR/wr.ssp370.2050.rds")sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp126 <- predict(sem.fennoscandia.5,wr.sf.126.2100,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp126 <- cbind(wr.sf.126.2100,wr.pred.ssp126)
wr.ssp126$TOC <- exp(wr.ssp126$fit)
wr.ssp126 <- merge(wr.ssp126, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp126,"WR/wr.ssp126.2100.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp370 <- predict(sem.fennoscandia.5,wr.sf.370.2100,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp370 <- cbind(wr.sf.370.2100,wr.pred.ssp370)
wr.ssp370$TOC <- exp(wr.ssp370$fit)
wr.ssp370 <- merge(wr.ssp370, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp370,"WR/wr.ssp370.2100.rds")The average TOC concentration was plotted.
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2050.rds")
temperature <- c(sem.pred.wr$Temp,wr.ssp126.2050$Temp,wr.ssp370.2050$Temp, wr.ssp126.2100$Temp,wr.ssp370.2100$Temp)
precipitation <- c(sem.pred.wr$Precip,wr.ssp126.2050$Precip,wr.ssp370.2050$Precip, wr.ssp126.2100$Precip,wr.ssp370.2100$Precip)
ndvi <- c(sem.pred.wr$NDVI,wr.ssp126.2050$NDVI,wr.ssp370.2050$NDVI, wr.ssp126.2100$NDVI,wr.ssp370.2100$NDVI)
runoff <- c(sem.pred.wr$Runoff,wr.ssp126.2050$Runoff,wr.ssp370.2050$Runoff, wr.ssp126.2100$Runoff,wr.ssp370.2100$Runoff)
tndep <- c(sem.pred.wr$TNdep,wr.ssp126.2050$TNdep,wr.ssp370.2050$TNdep, wr.ssp126.2100$TNdep,wr.ssp370.2100$TNdep)
toc <- c(sem.pred.wr$TOC95,wr.ssp126.2050$TOC,wr.ssp370.2050$TOC, wr.ssp126.2100$TOC,wr.ssp370.2100$TOC)
neworder <- c("1995","SSP126","SSP370")
model <- factor(c(rep("1995",dim(sem.pred.wr)[1]),rep("SSP126", dim(wr.ssp126)[1]), rep("SSP370", dim(wr.ssp370)[1])),levels = neworder)
compare.ssp <- data.frame(value = c(temperature, precipitation, ndvi, runoff,toc), variable = c(rep("Temperature (°C)", length(temperature)), rep("Precipitation (mm)", length(precipitation)), rep("NDVI",length(ndvi)), rep("Runoff (mm/y)",length(runoff)), rep("TOC (mg/L)", length(toc))), model = rep(model,5))
ggplot(compare.ssp)+geom_boxplot(aes(y= value,col=variable))+
labs(y="", x = "")+scale_color_brewer(type = "seq",palette = "Dark2",direction = -1)+
theme_minimal()+
theme(legend.position = "", axis.text.x = element_blank())+
facet_grid(cols = vars(model), rows = vars(factor(variable, levels = c("Temperature (°C)","Precipitation (mm)","NDVI","Runoff (mm/y)","TOC (mg/L)"))), scales = "free_y", shrink = T)
ggsave("WR/compare.ssp.png", width = 10, height = 20, units = "cm", dpi = "retina")
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2100.rds")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2100.rds")
list.sf <- list(wr.ssp126.2050, wr.ssp126.2100, wr.ssp370.2050, wr.ssp370.2100)
names(list.sf) <- c("SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
for(i in names(list.sf)){
list.sf[[i]]$diffTOC <- list.sf[[i]]$TOC - list.sf[[i]]$TOC95
list.sf[[i]]$diff.log <- list.sf[[i]]$fit - list.sf[[i]]$wr.pred.fit
list.sf[[i]]$exp.diff.log <- exp(list.sf[[i]]$diff.log)
}
all.diff <- data.frame(gid = c(wr.ssp126.2050$gid, wr.ssp370.2050$gid, wr.ssp126.2100$gid, wr.ssp370.2100$gid))
for(x in c("diffTOC","diff.log","exp.diff.log")){
all.vec <- c()
for(i in names(list.sf)){
vec <- select(st_drop_geometry(list.sf[[i]]),all_of(x))
all.vec <<- rbind(all.vec,vec)
}
all.diff <- cbind(all.diff, all.vec)
}
max.lim <- apply(all.diff, 2, max, na.rm = T)
min.lim <- apply(all.diff, 2, min, na.rm = T)
legends <- c("Difference TOC (mg/L)","Difference log(TOC)", "exp(difference log(TOC)")
diffs <- c("diffTOC","diff.log", "exp.diff.log")
g <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Difference log(TOC)\n SSP 1-2.6 - 2050", mid = 0, limits = c(min.lim[3],max.lim[3]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2050.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="b) Difference log(TOC)\n SSP 1-2.6 - 2100", mid = 0, limits = c(min.lim[3],max.lim[3]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2100.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="c) Difference log(TOC)\n SSP 3-7.0 - 2050", mid = 0, limits = c(min.lim[3],max.lim[3]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2050.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="d) Difference log(TOC)\n SSP 3-7.0 - 2100", mid = 0, limits = c(min.lim[3],max.lim[3]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2100.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")listgrob <- list("WR/SSP126.2050.diff.log.png", "WR/SSP126.2100.diff.log.png","WR/SSP370.2050.diff.log.png","WR/SSP370.2100.diff.log.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
ssp <- plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
save_plot("ssp.toc.png",ssp, ncol = 2, nrow = 2, base_height = 4, base_width = 3)
print(ssp)
4 Coastal drainage basins: coastal darkening
Most coastal drainage basins in Sweden and Finland drain into to the Baltic Sea (see Supplementary 2), while Norwegian sea watershed drain in Skagerrak, in the North sea, in the Norwegian Sea and in the Barent sea in the North (Sætre 2007). The drainage basin of the Baltic sea matches almost exactly the borders of Sweden and Finland : https://www.grida.no/resources/5324. All contribute to the Norwegian Coastal Current water masses. The estimated exported TOC was computed as:
\[ TOC_exp = [TOC] (mg.L^{-1}) \times Runoff (mm/y, or L.m^{-2}.y^{-1}) \times Watershed area (m^2)\]
Where the TOC concentration (past and future) was obtained by back-transforming the results of the model (log(TOC)) into TOC in mg/L. The result of the above equation, giving the amount of TOC exported in mg, is converted to Tg. Forecast are presented on Figure 15*.
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.centroid <- st_centroid(wr.sf.95, of_largest_polygon = T) %>% dplyr::select("gid","geometry")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
nor <- st_read("Country_shapefile/norway.shp") %>% st_transform(st_crs(wr.centroid))
swe <- st_read("Country_shapefile/sweden.shp") %>% st_transform(st_crs(wr.centroid))
fin <- st_read("Country_shapefile/finland.shp") %>% st_transform(st_crs(wr.centroid))
wr.norge <- wr.centroid[nor,]
wr.norge$Nation <- "Norway"
wr.sweden <- wr.centroid[swe,]
wr.sweden$Nation <- "Sweden"
wr.finland <- wr.centroid[fin,]
wr.finland$Nation <- "Finland"
wr.nation <- rbind(wr.norge,wr.sweden,wr.finland)
wr.sf.95 <- wr.sf.95 %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid") %>% sp::merge(dplyr::select(st_drop_geometry(sem.pred.wr), c("gid","wr.pred.fit","TOC95")), by = "gid")
names(wr.sf.95[which(names(wr.sf.95) == "TOC95")]) <- "TOC"
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")%>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2100.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2100.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
list.sf <- list(wr.sf.95, wr.ssp126.2050, wr.ssp126.2100, wr.ssp370.2050, wr.ssp370.2100)
names(list.sf) <- c("Data1995","SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
export <- function(sf){
discharge <- sf$Runoff*sf$area # Q in L/y
toc.export <- discharge * sf$TOC * 10^(-15) # amount of exported TOC in Tg/y
sf2 <- cbind(sf,discharge,toc.export)
return(sf2)
}
list.sf <- lapply(list.sf, export)
list.mean <- c()
for(i in 1:length(list.sf)){
list.mean[[i]] <- list.sf[[i]] %>% st_drop_geometry() %>% group_by(Nation) %>% summarise(total = sum(toc.export)) %>% as.data.frame()
}
names(list.mean) <- names(list.sf)
list.diff <- c()
for(i in 1:length(list.mean)){
df1 <- list.mean[[1]]
df2 <- list.mean[[i]]
diff <- df2[,2] - df1[,2]
x <- cbind(list.mean[[i]],diff)
names(x) <- c("Nation", names(list.mean)[i], paste("diff", names(list.mean)[i], sep = ""))
list.diff[[i]] <- x
}
toc.export.df <- list.diff[[1]] %>% merge(list.diff[[2]], by = "Nation") %>% merge(list.diff[[3]], by = "Nation") %>% merge(list.diff[[4]], by = "Nation") %>% merge(list.diff[[5]], by = "Nation") %>% melt(id.vars = c("Nation"), value.name = "TOC_Tg")
toc.export.df$Model <- NA
m.126 <- grep("126",toc.export.df$variable)
m.370 <- grep("370",toc.export.df$variable)
toc.export.df$Model[m.126] <- "SSP126"
toc.export.df$Model[m.370] <- "SSP370"
toc.export.df$Model[-c(m.126,m.370)] <- "Historical"
toc.export.df$Year <- NA
m.95 <- grep("95",toc.export.df$variable)
m.2050 <- grep("2050",toc.export.df$variable)
m.2100 <- grep("2100",toc.export.df$variable)
toc.export.df$Year[m.95] <- "1995"
toc.export.df$Year[m.2050] <- "2050"
toc.export.df$Year[m.2100] <- "2100"
toc.export.df$Diff <- NA
m.d <- grep("diff",toc.export.df$variable)
toc.export.df$Diff[m.d] <- "Difference"
toc.export.df$Diff[-m.d] <- "Absolute"
toc.export.sum <- toc.export.df %>% group_by(Model,Year,Diff) %>% summarise(TOC_Tg = sum(TOC_Tg))
toc.export.sum$Nation <- "Total"
toc.export.sum$variable <- NA
toc.export.df <- rbind(toc.export.sum,toc.export.df)
toc.export.df <- transform(toc.export.df, Nation = factor(Nation, levels = c("Norway","Sweden","Finland","Total")))
knitr::kable(toc.export.df) %>% kable_styling()
saveRDS(toc.export.df,"WR/toc.export.df.rds")
ggplot(filter(toc.export.df, Diff == "Difference")) + geom_point(aes(x=Year, y = TOC_Tg, col = Model, shape= Model)) +
geom_path(aes(x=Year, y = TOC_Tg, col = Model, group = Model), size = 0.3) +
labs(x="Year", y= "TOC (Tg)", subtitle = "b.")+
scale_color_manual(values = c("firebrick4","steelblue4","darkorange"))+
facet_grid(cols=vars(Nation))+
geom_abline(slope = 0, intercept = 0, col = "firebrick4", size = 0.3)+
theme_bw(base_size = 5)+
theme(legend.position = "bottom")
ggsave(filename = "WR/TOC_export_diff.png", dpi = "retina", width = 6, height = 6, units = "cm")
ggplot(filter(toc.export.df, Diff == "Absolute")) + geom_point(aes(x=Year, y = TOC_Tg, col = Model, shape= Model)) +
geom_path(aes(x=Year, y = TOC_Tg, col = Model, group = Model), size = 0.3) +
labs(x="Year", y= "TOC (Tg)", subtitle = "a.")+
scale_color_manual(values = c("firebrick4","steelblue4","darkorange"))+
facet_grid(cols=vars(Nation))+
geom_abline(slope = 0, intercept = 0, col = "firebrick4", size = 0.3)+
theme_bw(base_size = 5)+
theme(legend.position = "none")
ggplot(toc.export.df) + geom_point(aes(x=Year, y = TOC_Tg, col = Model, shape= Model)) +
geom_path(aes(x=Year, y = TOC_Tg, col = Model, group = Model), size = 0.3) +
labs(x="Year", y= "TOC (Tg)")+
scale_color_manual(values = c("firebrick4","steelblue4","darkorange"))+
facet_grid(cols=vars(Nation), rows = vars(Diff), scales = "free_y")+
geom_abline(slope = 0, intercept = 0, col = "firebrick4", size = 0.3)+
theme_bw(base_size = 8)+
theme(legend.position = "bottom")
ggsave(filename = "WR/TOC_export.png", dpi = "retina", width = 15, height = 10, units = "cm")
knitr::include_graphics("WR/TOC_export.png")
Figure 16: TOC exported from coastal drainage basins (mg/L)
References
CORDEX. 2021. “CORDEX experiment design for dynamical downscaling of CMIP6,” no. May: 1–8.
Detsch, Florian. 2021. “Download and Process GIMMS NDVI3g Data [R package gimms version 1.2.1].” Comprehensive R Archive Network (CRAN). https://cran.r-project.org/package=gimms.
“Environmental data NINAnor/NOFA Wiki.” n.d. Accessed February 21, 2022. https://github.com/NINAnor/NOFA/wiki/Environmental-data.
European Environment Agency. 2021. “Emissions of the main air pollutants in Europe INDICATOR ASSESSMENT.” https://www.eea.europa.eu/ims/emissions-of-the-main-air.
Grubb, Michael, Sujata Gupta, Kirsten Halsnaes, Bertjan Heij, Suzana Kahn Ribeiro, Shigeki Kobayashi, Mark Levine, and Daniel Martino. n.d. “Summary for Policymakers.” IPCC. https://doi.org/10.1017/9781009157896.001.
Hausfather, Zeke. 2019. “CMIP6: the next generation of climate models explained.” https://www.carbonbrief.org/cmip6-the-next-generation-of-climate-models-explained.
Hijmans, Robert J, Jacob Van Etten, Joe Cheng, Michael Sumner, Matteo Mattiuzzi, Jonathan A Greenberg, Oscar Perpinan, et al. 2018. “Package ’raster’ Type Package Title Geographic Data Analysis and Modeling.” https://github.com/rspatial/raster/issues/ https://cran.r-project.org/package=raster.
Myers-Smith, Isla H, Jeffrey T Kerby, Gareth K Phoenix, Jarle W Bjerke, Howard E Epstein, Jakob J Assmann, Christian John, et al. 2020. “Complexity revealed in the greening of the Arctic.” https://doi.org/10.1038/s41558-019-0688-1.
Raju, Komaragiri Srinivasa, and Dasika Nagesh Kumar. 2020. “Review of approaches for selection and ensembling of GCMS.” IWA Publishing. https://doi.org/10.2166/wcc.2020.128.
Sætre, Roald. 2007. The Norwegian Coastal Current - Oceanography and Climate. Edited by Tapir Academic Press. Trondheim.
The National Center for Atmospheric Research. 2018. “Global GIMMS NDVI3g v1 dataset (1981-2015).” http://poles.tpdc.ac.cn/en/data/9775f2b4-7370-4e5e-a537-3482c9a83d88/.
Tucker, Compton J, Jorge E Pinzon, Molly E Brown, Daniel A Slayback, Edwin W Pak, Robert Mahoney, Eric F Vermote, and Nazmi El Saleous. 2005. “An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data.” International Journal of Remote Sensing 26 (20): 4485–98. https://doi.org/10.1080/01431160500168686.